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Abstract 

We claim that the celebrated Stefan condition on the moving interphase, accepted in 
mathematical physics, can not be imposed if energy sources are spatially distributed in 
the volume. A method based on Tikhonov and Samarskii ideas for numerical solution of 
the problem is developed. Mathematical modelling of energy relaxation of some processes 
useful in modern ion beam technologies is fulfilled. Necessity of taking into account effects 
completely outside the Stefan formulation is demonstrated. 
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1 Introduction 



The Stefan problem concerns solid-liquid or liquid-vapor phase transitions when moving 
unknown beforehand surface S of phase transition is formed (see, e.g., In fact, a 

formulation of the Stefan problem was given for the first time by G. Lame and B.P. Clapeiron 
in 1831 for a particular case of equal temperature of liquid and crystalline phases |2J. In 
1889 J. Stefan published four papers devoted to the subject (in particular, to the description 
of soil freezing) in which the problem was formulated in a general form According to it, 
for the interphase the following condition 

dT(x s + 0,t) dT(x s -0,t) 
ol ~dx q ~dx = L PsoiVs, (1) 

defining Stefan's problem, has been suggested. Here Vs — d£ s /dt is the velocity of the 
boundary surface S, K so i and Ku q are thermal conductivities of material for solid and liquid 
phases, L and p so i are the melting heat and density, correspondingly. Condition (1) has a 
clear physical meaning. Indeed, according to the Fourier law, heat flow j is proportional to 
the temperature gradient, 

j = —K grad T. 

Therefore, the left-hand side of (1) is the heat absorbed in the unit of area per the unit of 
time. The expression in the right-hand side is the heat connected with freezing or melting 
of material crossed per the unit of time by the the unit of area. 

A complete mathematical formulation of the Stefan problem includes, besides (1), a 
condition of continuity on the surface S separating solid and liquid phases 

T| S = T*, (2) 

where T* denotes temperature of the phase transition, and the energy conservation law 

dT 

pC — = -div j + q(x,t). 

Here q(x, t) represents the power of external heat sources, C is the specific heat. In the 
original Stefan papers q(x,i) = 0, so that the whole heat transfer has been considered to be 
a consequence of the temperature gradient inside the medium. 

If one also specifies initial and boundary conditions, the Stefan problem can be solved 
more often approximately, but sometimes exactly. Particular examples of suitable boundary 
conditions are considered below. 

Relations (1) and (2) are usually used in numerical algorithms explicitly. Another ap- 
proach was suggested by A.N. Tikhonov and A. A. Samarskii in 1953 @|. According to it, 
conditions (1) and (2) are included themselves into the energy conservation equation to 
obtain generalized formulation of the Stefan problem in the form 

(pC + LS{T- T*)) + v grad T^j = div(K grad T) + q(x, t), (3) 

where the term L S(T — T*) dT/dt describes the additional heat input expended on the 
phase transformation, v grad T takes into account possible temperature change due to 
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convection (hereafter we ignore it for simplicity). The main idea of this approach is quite 
clear, too. Namely, it is suggested to treat the heat of fusion L as an additional component 
of the thermal capacity pC which gives contribution only at the point of phase transition. 

Lately Samarskii and his followers have turned this idea into effective numerical al- 
gorithms (see, e.g. But even in those papers equation (3) is considered only as a 
corollary of the condition (1). For example, it was derived in pp by substituting expression 
L S(T — T*) dT/dt instead of the term L S(x — £s(i)) Vs, which is assumed to be included 
in the heat equation to account for the heat absorption on 2-dimensional interface S. 

The purpose of this paper is to show that the condition (3) supplies us with a more 
powerful description of phase transitions, that may be used even in the case when (1) and 
(2) are not applicable. 



2 Heuristic arguments 

As it was mention above, the possibility of solving the classical Stefan problem by making 
use of condition (3) has been demonstrated by Samarskii and his co-authors. Therefore, we 
only consider an example when (3) is applicable and (1), (2) are not. To this end let us 
study the following problem: 

dT 

(pC + LS(T- T*)) -fodiv{k grad T) + q(t), (4) 
T(x,0) =T <T*, 

where all parameters of (4) are suggested to be independent of x. Due to the spatial 
uniformity, it is evident that the condition 

grad T = 

holds on the solutions of (4). In this case Eq. (4) is reduced to an ordinary differential one 

(pC + LS(T-T*))^ = q(t) (5) 

with the initial condition 

T(0) = T . 

Integrating both sides of (5) over t just near the phase transition temperature T*, one 
obtains 

/ (pC + LS{T ~T*))dT = / q(t)dt, (6) 
Jt'-o Jt 

where St is a time necessary for the phase transition. It is evident from (6) that 

St > (7) 

where Q is the maximum value of q(t) in the interval (t,t + St). The inequality (7) means 
that the phase transition at a fixed spatial point lasts a finite, distinct from zero, time. 

This simple example shows something completely different from the Stefan description 
of the phase transition. Let us examine it carefully. 
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1. First of all, instead of gradual warming (or cooling) up the pattern due to the influence 
of one of its boundaries, here we have an uniformly heated layer. Therefore, creation of 
2-D surface S(y, z) separating the solid and liquid phases in x is evidently impossible 
due to a total equivalence of all spatial points x. 

2. One can also expect that finiteness of the phase transition time, St, forces all points 
within a spatial layer of nonzero thickness to be at the same temperature T*. 
This is expected even in the case when the power deposition q(x,t), unlike in the 
example considered, is spatially irregular 1 . 



T(x) 





a) b) 

Fig. 1. Evolution of the temperature distribution in the pattern for: a) a case 
of spatially distributed sources of heat, b) the classical Stefan problem (with 
heating from left to right). Here t\ < t 2 < t 3 . 

If we consider (^-function in (3) as a limit of a bounded function D(T — T*) localized in 
the vicinity of T — T*, then the possibility to obtain the solution shown in Fig. 1(a), or 
analytically 

i)T 

0, grad T — ► 0, 



dT 



follows from the indefinitencss 



D{T -T*) 



dT 

~dt 



oo 



0. 



springing up in the left-hand side of (3). Clearly, that this indcfinitcness can take a finite 
value and compensate in a space region with nonzero thickness the spatially distributed 
source q(x,t) which contributes to the right-hand side of (3). This, of course, is no more 
true if the external sources are absent and heat enters the pattern only through its boundary. 

1 Indeed, let material has just reached the temperature T* at some point x and now begins receiving its 
portion of the heat necessary for melting. Then another adjacent point x + Ax, which attained the melting 
temperature merely a little earlier, can be still in the state of heat receiving and, therefore, must have the 
same temperature T* (see Fig.l). 
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In a general case, one can expect existence of two jumps for spatial derivatives of 
temperature on the boundaries S of the volume Vt* with T = T*, instead of one for the 
classical Stefan problem, but the condition (1) is hardly met for any of them (see Fig. 1(a), 
where intersection of the boundary S by (x,y)-plane in points 1, 2, 3 and 4 is seen). Indeed, 
to prove the existence of two jumps — one from the side of the solid and another from the 
side of the melted phase — it is sufficient only to show that the spatial derivative on the 
surface S, taken externally, is not equal to zero. The co-ordinates of the boundary £(t) 
can be found as a solution of an equation 

T(x, t) - T* = 0, 

where T(x, t) is the solution of the heat equation (3) outside the volume Vt- ■ Taking the 
total temporal derivative, one obtains 

— +gradT ■ - = 0. 

Thus grad T = automatically implies dT/dt = 0. It is evident that such conditions are 
impossible if the external sources are not adjusted specially to stabilize the temperature 
in the infinitesimal layers adjacent to the volume Vt* just before and just after the phase 
transition. 



3 Beam induced phase transitions 

To verify the conclusions which we have just come to, let us study numerically the dynamics 
of phase transition induced by a short powerful ion beam in solids. At present this technology 
is really used for modification of surface layers to create new materials with unique physical 
and chemical properties (see, e.g. The process is underlain by the equation for heat 

transfer which we discussed in the previous sections: 

dT d ( dT\ 
PiT)c{T)°4 = 4- (k{T)°— )+ q . (8) 



dt dx \ dx r 
The initial and boundary conditions could be taken in the form: 

rpl n X rp dT(0,t) dT(h,t) 

Let us consider, for definitcness, an iron pattern, which thermal properties are described 
in popular reference books, and choose dimcnsionlcss (DL for brevity) variables 

T := T/Tq, x :— x/Iq, t := t/r 

as follows: 

T = 293 K, l a = 1CT 5 m (the pattern thickness), 

r = 3 10~ 7 s (duration of ion beam pulse from an accelerator). 

For DL power deposition q we take a simple model, shown in Fig. 2 and 3, 
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q(x = 0,i) q{x,t = 0.5) 




t x 
Fig. 2. Power deposition: i-dependence Fig. 3. Power deposition: x-dependence 



with analytical representation 

q(x,t) = Q qi(x)q 2 (t), 

where 

n{ z ) = ^j— 7 \ 

1 + exp^j(z - Zi) 

and Q describe the total DL energy brought into the pattern (here Q — 59.44, x\ = 
0.07, t\ = 1, fii = 100). For simplicity, we neglect in (8) a small difference between physical 
parameters for the solid and liquid phases. 

Now, using of the general idea due to Tikhonov and Samarskii 0], we assume an expres- 
sion: 

p{T)c(T) = 1 + X5(T - T*, A) 

for DL specific heat, where A denotes the DL heat of fusion and 5(T — T*, A) is an approx- 
imate 5-function, smoothed with the help of the Gaussian distribution of width A (see Fig. 
4) 2 . 

Now Eq. (8) can be solved numerically on the space-time grid x and t with steps h x and 
h tl within intervals x £ (0, 1), t £ {0,t max ): 

Xj =h x -j, j = 0, . . . , n x , h x = l/n x , 
tk = ht-k, k = 0, . . . , n t , h t — t max /n t , 
where n x and n t are numbers of partitions. 

2 There were other methods of smoothing in original papers by Samarskii et al. They used regularization 
on the space grid. 



G 




T 



Fig. 4. Approximate <5-function with maximum 
at T* and smearing equals A (for A = 0, 03) 

The following difference scheme with weights 7 was implemented (see [7] for details): 



T 

„k 3 



fe+1 



rpk+1 



2T fc + l 



rpk + 1 

1 3-l 



hi 



+ (l- 7 ) 



hi 



(9) 



q(xj,t k 



— ) 

2 J ' 



where 

Tf = T( Xj ,t k ), e) = p{Tf)c{Tf), 

and the upper index numerates different moments of time (time "levels"), the lower one 
specifies a set of spatial co-ordinates. The scheme is absolutely convergent at 7 = 0.5 and 
possesses the second-order accuracy for both variables. 

From initial condition T(x, 0) = To, values T® (J = 0, . . . , n x ) on a zero time level are 
known. The boundary conditions 



rpk 



rpk 



= 0, k = 1, 



, n, 



allow one to introduce symmetric points X-\ — —h x and x nx +\ = 1 + h x with appropriate 
values T*! = T* and +1 = T% _ x respectively. So, we can use the Eq. (9) in points Xo 
and a;,^ . Using initial and boundary conditions, we obtain a system of n x linear algebraical 
equations with the same number of variables. Thus, under the accepted approximation, we 
reduced the partial differential equation (8) to system (9) of linear algebraic equations. The 
matrix of this system is tridiagonal and after its solution 3 we obtain value Tj (J = 0, . . . , n x ) 
at the first time level. Repeating this process, values on every time level k are computed. 

The result of straightforward verification of the Stefan condition (f ) is shown in Fig. 5, 
where the function 

dT 

dx 

v + A + ti / 

3 Recursive relations for determining the solution of algebraic problem (9) comprise the well-known sweep 
method, called also forward-backward or Thomas algorithm 1 . 



dx 



4 

dt 



(10) 
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is depicted. 



0.1 0.2 0.3 0.4 0.5 

t 

Fig. 5. 



The left and right points, in which the spatial derivatives of temperature were taken 
in (10), are shown in Fig. 4. They define a spatial layer which nearly the whole fusion 
energy is absorbed within. From Fig. 5 one can see that condition (1) is satisfied indeed, 
but only after a characteristic relaxation time tx has elapsed. The physical meaning 
of t\ is clear from Fig. 6. Namely, it corresponds to the transition from a rapid to slow 
motion of the exterior interphase surface. In the case when boundary motion is rapid, the 
heat necessary for fusion is brought into the melting layer directly from the external 
source q(x, t). The slow motion corresponds to the ordinary Stefan mode when the process 
is controlled mainly by the heat entered into the layer through its boundary. 
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It is also seen from Fig. 6 that transition to the Stefan mode takes place earlier than the 
external source to be totally turned off: 



h < t. 

Time ti shown in Fig. 6 denotes a moment when the thickness of the melted material begins 
to diminish due to heat escape into a more cooler solid phase. 

Fig. 7 and 8 also confirm the conclusions which we have come to in the previous section. 
Formation of the "tableland" (whose height corresponds to the fusion temperature) for 
spatial temperature distribution is distinctly seen in Fig. 7. The narrow strip restricted 
by two dashed lines in Figs. 7 and 8 exhibits the width of the smoothed <5-function. We 
believe that existence of two breaks for the spatial derivative is masked in Fig. 7 with this 
(5-function smearing. Fig. 8 demonstrates a temperature evolution for two divorced spatial 
points. One can make sure that the above mentioned time interval corresponding to the 
same temperature at the different spatial points really exists. It is evident that such a 
behavior of temperature has nothing to do with the traditional description in the framework 
of (1) and (2). 





Fig. 7. Spatial temperature distribution 
for: 1) t = 0,162; 2) t = 0,186; 3) t = 
0,204; 4) t = 0,216 (A = 0,01) 



Fig. 8. Evolution of temperature for: 
1) x = 0; 2) x = 0,04 (A = 0,01) 



Fig. 9 shows a time-dependence of the interphase coordinate. Numbers 1 and 2 denote 
the regions where verification of the Stefan condition (1) is impossible due to A-instability. 
It means that small variations of fusion temperature value, T* , lead to a drastic change of 
interphase position (see dotted lines in Fig. 9). 
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1 



2 



3 



4 



t 



Fig. 9. Time dependence of the interphase coordinate. Here a 
curve in the middle corresponds to the fusion temperature equals 
T* . For the upper and lower curves, the fusion temperature is 
chosen to be at T* - A/2 and T* + A/2, accordingly (A = 0, 01) 



4 Track formation in solids 

The next example demonstrating the preference for the (^-function approach is connected 
with the problem of track formation in solids. In fact, at present nobody knows with 
certainty the main mechanism responsible for these processes. Furthermore, it seems like 
the universal model explaining all of them does not exist and different materials behave 
differently under heavy ion attack. Here we assume the so-called thermal spike model based 
on the following system of two coupled nonlinear differential equations (see, e.g. |Hj and 
references therein): 



where T e and Ti are electrons and lattice temperatures, respectively, C e , Cj and K e ,Ki 
specific heat and thermal conductivity for the electronic system and lattice, p is the material 



P C e (T e )-^ = -- rK e (T e )-^ -fl-(T e -T0 + «(r,i), 



(11) 




(12) 
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density, g the electron-atom coupling, q(r, t) the power brought on the electronic system, r 
the radius in cylindrical geometry with the ion path as the axis. One can see that electrons 
receive their energy directly from the external source q(r, t) which takes into account ion 
energy loss in electron gas. The characteristic duration of source activity is usually in the 
range 1CP 15 — 5 x 10~ 15 s. According to (12), atoms are heated due to electron-atom coupling 
represented by the term g ■ (T e — Ti). Nuclear interaction of atoms with the projectile ion 
is relatively small and, therefore, can be neglected. It is clear that coupling is the most 
effective at the beginning of the relaxation process when T e ^> Tj and g ■ (T e — T,) ~ gT e . 
The initial conditions can be chosen in a form 

T e (r,0) =Ti(r,0) =T , 

and the boundary ones 4 can be taken as 

r=0 V 9r J r=0 

where r max was taken of order 10 -5 cm. 

The thermal spike model explains track formation as a structural transition of lattice 
due to its warming-up and melting with subsequent quenching. These processes are usually 
accompanied with disorder creation in the lattice. Indeed, rapid quenching leads to a "con- 
servation" of atoms' random places that were in the melted material just before cooling. 
For amorphous materials, which are characterized by high disorder of atoms' positions and 
small values of thermal conductivity, quenching, quite the contrary, leads to putting atoms' 
places in order. But in either case, structural modifications are observed in the microscope 
as an ion trace in solid. 

Besides thermal spike, one may assume the ion spike as well, when the track is formed 
due to the electrostatic repulsion of ionized atoms. The main reason justifying our utilization 
of system (11), (12) is an agreement of nuclear track radii, calculated in this framework, 
with the experimental data [5| . The total formulation of the model includes many physical 
details, such as a description of the source q(x, t), and is outside the scope of this publication. 
Here we only touch some problems concerning the main topic of the paper. 

A numerical algorithm similar to that described above has been elaborated for numerical 
solving system (11), (12). The radial distribution of the lattice temperature Ti around the 
path of Pb in amorphous Ge at kinetic energy of impinging ions of about 110 MeV is shown 
in Fig. 10 (for two different moments of time). One can see the typical "tablelands" similar 
to those discussed in the previous chapter and which could not be obtained in frames of 
the classical Stefan approach. Our calculations show that the "tableland" exists here only 
during a transitory time t\, when the material is under a strong exposure of the source 
g ■ (T e — Ti) ~ g T e . It is shorter than r = gC e /g f» 10~ 12 s (see Fig. 10), where r is 
electron-atom relaxation time 5 , characterizing duration of source activity in (12) (compare 
this conclusion with data presented in Fig. 6). 

4 One should take into account that there is no heat transfer at the center of track. 
5 The formula for r estimation follows from Eq. (11). 




11 



T[°C] 

1800r 




2 4 6 8 10 12 14 16 18 20 



R [nm] 

Fig. 10. The radial distribution of the lattice temperature Ti 
for two different moments of time 

It is interesting to note that there is a real, met in the nature, "regularization" of 
5-function analogous to that implemented in this paper. For materials with a complex 
molecular structure (high temperature superconductors, biological molecules, alloys etc.), 
the melting temperature is not fixed but, instead, smeared within a characteristic interval 
where atom bonds of different type are gradually destroyed with temperature increase. In 
this case the only possible approach to the problem should be based on the condition (3). 
An approximate ^-function, analogous to that shown in Fig. 4, can be extracted here directly 
from the experiment. In 9 , a model based on the smeared 5-function approach and Eqs. 
(11), (12) was used for computation of effective electron-atom relaxation time r in a high 
temperature superconductor. The established t turned out to be in a good agrement with 
experimentally observed values. 

5 Conclusion 

To the best of our knowledge, the peculiarities of phase transition dynamics, we discussed 
in this paper, have never been considered explicitly in mathematical physics. This fact 
may be explained partially by the necessity to use very powerful spatially distributed 
external sources of heat, in order the above mentioned effects to be urgent. Such sources 
were hardly available for industrial applications not long ago. However, the examples which 
given above are likely evidences of the fact that such sources, "interfering" in the thermal 
conductive processes, are integral parts of all most recent ion beam technologies. The 
numerical investigations, which have been undertaken, show that the (5-function approach 
to phase transitions is a suitable instrument to tackle these problems, though the authors 
of this idea have never used it in such a context. 

Within networks of scientific papers devoted to the interphase motion problem, one can 
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distinguish in retrospect the following logical order: formulation of the Stefan problem (Lame 
and Clapeiron; Stefan) — ► application of the ^-function approach for equivalent representa- 
tion of it (Tikhonov and Samarskii) — ► numerical implementations of the idea (Samarskii 
with co-authors; this review) — ► description of materials with interval distributed fusion 
temperature (a natural physical interpretation of the previous step). Here we have shown 
that it is more expedient to turn over this order and take its last element as the basis for 
solving both the classical Stefan problem and a more general one including spatially dis- 
tributed sources. In other words, both of these solutions could be considered as an idealized 
limiting case A — > of a natural physical point of view that none of phase transitions take 
place at the exactly defined value of fusion temperature. A peculiarity of a new, found in 
this paper, solution to Eq. (3) is its "tableland" behavior seen in Figs. 1(a), 7, 8 and 10. 
Such solutions could never be obtained in the framework of the classical Stefan formulation 
(compare with Fig. 1(b)). 
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